An Reproducible Example of Colocalization Analysis
Yonghui Dong
10/19/2020
(1) Load packages and read data
The dataset is downloaded from Metabolights. Accession number: MTBLS487. The data was obtained in positive ion mode over an m/z range of 300–2000 using a Thermo MALDI LTQ Orbitrap XL instrument (Thermo Fisher Scientific, Bremen, Germany). The spatial resolution was 50 micrometers. For more information about the dataset, readers are kindly referred to the original publication (Hall et al., 2017).
library(Cardinal)
library(BiocParallel)
library(plotly)
path <- paste0("Data/test_POS", ".imzML")
Brain <- readMSIData(path, resolution = 5, units = "ppm",
mass.range = c(300, 2000), attach.only = T)(2) Pre-processing
Brain2 <-
Brain %>%
normalize(method = "tic") %>%
peakPick(method = "mad", SNR = 6) %>%
peakAlign(tolerance = 5, units = "ppm") %>%
process(BPPARAM = SerialParam())Check the ion image of m/z 844.527 before and after pre-processing in order to make sure that chosen pre-processing parameters are fine.
darkmode()
##(2.1) before pre-processing
image(Brain, mz = 844.527, smooth.image = "gaussian", plusminus = 0.03, colorscale = magma,
contrast.enhance = "suppression", normalize.image = "linear")
##(2.2) after pre-processing
image(Brain2, mz = 844.527, smooth.image = "gaussian", plusminus = 0.03, colorscale = magma,
contrast.enhance = "suppression", normalize.image = "linear")
(3) Co-localization
Here first 100 colocalized ions are retained.
#(3.1) co-localization PC36:1
coloc_826 <- colocalized(Brain2, mz= 826.572, n = 100, BPPARAM = SerialParam())
#(3.2) co-localization PC38:6
coloc_844 <- colocalized(Brain2, mz= 844.525, n = 100, BPPARAM = SerialParam())
#(3.3) co-localization PC40:6
coloc_872 <- colocalized(Brain2, mz= 872.557, n = 100, BPPARAM = SerialParam())(4) Extract in-source fragments, and create pseudo MS/MS spectrum
#(4.1) PC36:1
mycoloc <- as.data.frame(coloc_826)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw() +
ggtitle("m/z 826.577")
ggplotly(p)#(4.2) PC38:6
mycoloc <- as.data.frame(coloc_844)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw() +
ggtitle("m/z 844.528")
ggplotly(p)#(4.3) PC40:6
mycoloc <- as.data.frame(coloc_872)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw()+
ggtitle("m/z 872.559")
ggplotly(p)(5) Plot MALDI images
#(5.1) PC36:1
mycoloc <- as.data.frame(coloc_826)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_826.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()
#(5.2) PC38:6
mycoloc <- as.data.frame(coloc_844)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_844.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()
#(5.3) PC38:6
mycoloc <- as.data.frame(coloc_872)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_872.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()Session information
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.2 (2020-06-22)
## os macOS Catalina 10.15.7
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Asia/Jerusalem
## date 2020-10-19
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.0.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.0)
## biglm 0.9-2 2020-05-26 [1] CRAN (R 4.0.0)
## Biobase 2.48.0 2020-04-27 [1] Bioconductor
## BiocGenerics * 0.34.0 2020-04-27 [1] Bioconductor
## BiocParallel * 1.22.0 2020-04-27 [1] Bioconductor
## bitops 1.0-6 2013-08-17 [1] CRAN (R 4.0.0)
## Cardinal * 2.6.0 2020-04-27 [1] Bioconductor
## cli 2.1.0 2020-10-12 [1] CRAN (R 4.0.2)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 4.0.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.0)
## crosstalk 1.1.0.1 2020-03-13 [1] CRAN (R 4.0.0)
## data.table 1.13.0 2020-07-24 [1] CRAN (R 4.0.2)
## DBI 1.1.0 2019-12-15 [1] CRAN (R 4.0.0)
## digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.0)
## dplyr 1.0.2 2020-08-18 [1] CRAN (R 4.0.2)
## EBImage * 4.30.0 2020-04-27 [1] Bioconductor
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.0)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.0)
## farver 2.0.3 2020-01-16 [1] CRAN (R 4.0.0)
## fftwtools 0.9-9 2020-10-03 [1] CRAN (R 4.0.2)
## generics 0.0.2 2018-11-29 [1] CRAN (R 4.0.0)
## ggplot2 * 3.3.2 2020-06-19 [1] CRAN (R 4.0.0)
## glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.0)
## htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.0)
## htmlwidgets 1.5.2 2020-10-03 [1] CRAN (R 4.0.2)
## httr 1.4.2 2020-07-20 [1] CRAN (R 4.0.2)
## irlba 2.3.3 2019-02-05 [1] CRAN (R 4.0.0)
## jpeg 0.1-8.1 2019-10-24 [1] CRAN (R 4.0.0)
## jsonlite 1.7.1 2020-09-07 [1] CRAN (R 4.0.2)
## knitr 1.30 2020-09-22 [1] CRAN (R 4.0.2)
## labeling 0.3 2014-08-23 [1] CRAN (R 4.0.0)
## lattice 0.20-41 2020-04-02 [1] CRAN (R 4.0.2)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.0.0)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.0)
## locfit 1.5-9.4 2020-03-25 [1] CRAN (R 4.0.0)
## magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.0)
## MASS 7.3-53 2020-09-09 [1] CRAN (R 4.0.2)
## Matrix 1.2-18 2019-11-27 [1] CRAN (R 4.0.2)
## matter 1.14.0 2020-04-27 [1] Bioconductor
## mclust 5.4.6 2020-04-11 [1] CRAN (R 4.0.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.0)
## nlme 3.1-149 2020-08-23 [1] CRAN (R 4.0.2)
## pillar 1.4.6 2020-07-10 [1] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.0)
## plotly * 4.9.2.1 2020-04-04 [1] CRAN (R 4.0.2)
## png 0.1-7 2013-12-03 [1] CRAN (R 4.0.0)
## ProtGenerics * 1.20.0 2020-04-27 [1] Bioconductor
## purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.0)
## R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.0)
## Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2)
## RCurl 1.98-1.2 2020-04-18 [1] CRAN (R 4.0.0)
## rlang 0.4.8 2020-10-08 [1] CRAN (R 4.0.2)
## rmarkdown 2.4 2020-09-30 [1] CRAN (R 4.0.2)
## S4Vectors * 0.26.1 2020-05-16 [1] Bioconductor
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.0)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
## signal 0.7-6 2015-07-30 [1] CRAN (R 4.0.0)
## sp 1.4-4 2020-10-07 [1] CRAN (R 4.0.2)
## stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.0)
## tibble 3.0.3 2020-07-10 [1] CRAN (R 4.0.2)
## tidyr 1.1.2 2020-08-27 [1] CRAN (R 4.0.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.0)
## tiff 0.1-5 2013-09-04 [1] CRAN (R 4.0.0)
## vctrs 0.3.4 2020-08-29 [1] CRAN (R 4.0.2)
## viridisLite 0.3.0 2018-02-01 [1] CRAN (R 4.0.0)
## withr 2.3.0 2020-09-22 [1] CRAN (R 4.0.2)
## xfun 0.18 2020-09-29 [1] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library